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Abstract 

High  order  path-conservative  schemes  have  been  developed  for  solving  nonconservative 
hyperbolic  systems  in  [32,  7,  6].  Recently,  it  has  been  observed  in  [3]  that  this  approach  may 
have  some  computational  issues  and  shortcomings.  In  this  paper,  a  modification  to  the  high 
order  path-conservative  scheme  in  [7],  based  on  the  high  order  finite  volume  WENO  scheme 
with  subcell  resolution  and  utilizing  the  exact  Riemann  solver  to  catch  the  right  paths  at  the 
discontinuities,  is  proposed  to  improve  its  computational  performance  and  to  overcome  some 
of  the  shortcomings.  An  application  to  one- dimensional  compressible  two-medium  flows  of 
nonconservative  or  primitive  Euler  equations  is  carried  out  to  show  the  effectiveness  of  this 
new  approach. 
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1  Introduction 


Recent  years  have  seen  a  growing  interest  in  developing  numerical  algorithms  for  solving 
compressible  multicomponent  flows.  The  dynamics  of  inviscid  multicomponent  fluid  may 
be  modeled  by  the  Euler  equations.  However,  computations  often  run  into  unexpected 
difficulties  due  to  nonphysical  oscillations  generated  at  the  vicinity  of  the  material  interface, 
while  such  oscillations  do  not  arise  in  single-fluid  computations  when  a  nonlincarly  stable 
scheme,  such  as  the  essentially  non-oscillatory  (ENO)  or  weighted  ENO  (WENO)  scheme 
[16,  38,  29,  20],  is  used.  The  underlying  mechanisms  have  been  analyzed  and  several  methods 
have  been  developed  to  overcome  these  difficulties,  e.g.  in  [23,  22,  10,  1,  2,  25,  5]. 

There  are  mainly  two  approaches  to  circumvent  these  oscillations.  One  is  still  based  on  the 
conservative  Euler  equations,  and  the  other  is  to  write  the  Euler  equations  in  nonconservative 
or  primitive  form. 

The  ghost  fluid  method  (GFM)  developed  in  [13]  with  the  isobaric  fix  technique  in  [14] 
has  provided  an  attractive  and  flexible  way  to  treat  the  two-medium  flow  model  for  conser¬ 
vative  Euler  equations.  The  GFM  treats  the  material  interface  as  an  internal  boundary,  and 
by  defining  ghost  cells  and  ghost  fluids,  the  two-medium  flow  can  be  solved  via  two  respec¬ 
tive  single-medium  GFM  Riemann  problems.  This  method  is  simple  and  it  easily  extends 
to  multi-dimensions,  and  it  can  maintain  a  sharp  interface  without  oscillations.  Variants 
of  the  original  GFM  and  their  applications  can  be  found  in,  e.g.  [25,  26,  27,  46]  and  the 
references  therein.  Later,  these  techniques  are  used  to  develop  the  Runge-Kutta  discontinu¬ 
ous  Galerkin  (RKDG)  finite  element  method  for  multi-medium  flow  in  [33,  34,  51,  47].  The 
GFM,  even  though  approximating  directly  the  conservative  Euler  equations,  is  in  general 
not  a  conservative  method  because  the  flux  at  the  interface  is  double-valued.  Recently,  a 
conservative  modification  to  the  GFM  using  the  fifth  order  finite  difference  WENO  scheme 
with  third  order  Runge-Kutta  time  discretization  has  been  studied  in  [28] ,  which  attempts 
to  reduce  the  conservation  error  of  the  GFM  without  affecting  its  performance. 

The  other  approach  is  based  on  the  observation  that  erroneous  pressure  fluctuations 
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are  generated  by  the  conservative  equations,  hence  a  better  approximation  can  be  obtained 
when  writing  the  equations  in  a  nonconservative  (primitive)  form  [21,  22],  In  these  papers 
a  scheme  for  the  nonconservative  Euler  formulation  is  proposed,  with  consistent  correction 
terms  to  remove  the  leading  order  conservation  errors.  This  scheme  can  completely  eliminate 
spurious  oscillations  at  the  material  interface,  yielding  clean  monotonic  solution  profiles. 
This  work  has  been  extended  to  solve  two  dimensional  problems  in  [35].  In  this  method, 
the  correction  terms  depend  heavily  on  the  corresponding  conservative  numerical  scheme, 
and  only  second  order  nonconservative  scheme  has  been  provided  and  applied  to  shocks  of 
weak  to  moderate  strengths,  it  is  less  justifiable  for  high-resolution  schemes  with  narrow 
shock  transition  and  it  is  not  justified  in  cases  of  strong  shocks  [18].  It  is  noted  in  these 
papers  that,  for  nonconservative  hyperbolic  systems,  the  shock  relationships  are  not  uniquely 
defined  by  the  limiting  left  and  right  states,  but  also  depend  on  the  viscous  path  connecting 
the  two  states.  The  correct  shock  capturing  lies  in  getting  correctly  the  viscous  path. 

The  theory  developed  by  Dal  Maso,  LeFloch  and  Murat  [11]  gives  a  rigorous  definition  of 
nonconservative  products,  associated  with  the  choice  of  a  family  of  paths.  Later,  people  have 
paid  much  attention  to  the  development  of  numerical  schemes  for  solving  nonconservative 
hyperbolic  systems,  see  for  example  [32,  7,  8,  6]  and  references  therein.  A  high  order  Roe-type 
scheme  based  on  the  reconstructed  states  has  been  provided  in  [7]  for  the  one-dimensional 
case,  and  then  extended  to  two  dimensions  in  [6],  but  only  with  applications  to  shallow- 
water  systems.  Other  work  based  on  this  high  order  method  has  been  carried  out  for  solving 
two-phase  flow  models  in  [45,  12,  42],  A  limitation  of  this  approach  has  been  pointed  out 
recently  in  [3]  when  applied  to  nonconservative  Euler  equations.  The  problem  is  related  to 
the  effective  choice  of  a  correct  path  in  the  nonconservative  high  order  scheme.  It  appears 
that,  for  nonconservative  Euler  equations,  using  the  Roc  linearization  to  choose  a  path  as  in 
[43,  7,  6]  might  end  up  in  converging  to  a  weak  solution  of  a  different  path  [3]. 

In  this  paper,  we  focus  on  adapting  high  order  schemes  for  solving  nonconservative  Euler 
equations,  still  based  on  high  order  Roe- type  finite  volume  schemes  as  those  in  [7,  6].  For 
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the  nonconservative  or  primitive  Euler  equations,  the  right  path  should  recover  the  weak 
solution  of  the  conservative  Euler  formulation  with  density,  momentum  and  energy  as  its 
variables.  In  smooth  regions,  the  conservative  and  nonconservative  Euler  equations  are 
equivalent,  but  at  discontinuities,  they  are  not  [21].  As  is  well  known,  shock  capturing 
schemes  such  as  monotone,  total  variation  diminishing  (TVD),  or  ENO  and  WENO  schemes, 
smear  discontinuities  with  one  or  several  transition  points.  These  transition  points  are 
necessary  for  conservation,  however  in  a  nonconservative  formulation,  they  may  not  land  on 
the  correct  path  and  hence  may  lead  to  convergence  to  erroneous  weak  solutions  on  different 
paths.  Realizing  this  difficulty,  which  unfortunately  is  generic  with  all  shock  capturing 
schemes,  our  basic  idea  in  this  paper  is  to  use  Harten’s  subcell  resolution  technique  [15] 
to  sharpen  the  discontinuities  and  effectively  eliminate  (or  at  least  significantly  reduce)  the 
transition  points.  As  a  result,  the  convergence  towards  the  correct  weak  solution  based  on 
the  originally  desired  path  seems  to  be  restored. 

Harten’s  subcell  resolution  idea  [15]  is  based  on  ENO  schemes  with  a  Lax-Wendroff  time 
discretization  procedure  in  a  cell-averaged  framework.  Later,  this  idea  is  extended  in  [39]  to 
both  finite  volume  and  finite  difference  ENO  schemes  with  Runge-Kutta  time  discretization. 
Recently,  this  subcell  resolution  idea  has  been  used  in  solving  advection  equations  with  stiff 
source  terms,  to  obtain  correct  shock  speed  on  coarse  meshes  [48].  In  this  paper,  with  the 
sharp  left  and  right  states  at  the  discontinuities,  we  use  the  exact  Riemann  solution  to  catch 
the  right  path  that  connects  the  two  states,  as  the  correct  capturing  of  the  shock  speed  is 
sensitive  to  the  accuracy  of  the  numerically  achieved  path. 

In  Section  2,  we  first  describe  the  high  order  Roe  scheme  for  the  nonconservative  hy¬ 
perbolic  system.  Then  we  introduce  the  two  medium  flow  model  for  nonconservative  Euler 
equations  in  Section  3.  In  Section  4,  we  illustrate  how  to  apply  the  WENO  reconstruction 
with  subcell  resolution  to  the  high  order  Roe  scheme.  In  Section  5,  we  describe  the  compu¬ 
tation  of  the  integral  term  in  the  high  order  Roe  scheme  in  detail.  In  Section  6,  we  present 
the  level  set  function  to  track  the  material  interface.  A  summary  of  our  algorithm  to  aid 
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implementation  is  given  in  Section  7,  and  one-dimensional  numerical  examples  are  provided 
to  demonstrate  the  effectiveness  of  our  approach  in  Section  8.  Concluding  remarks  follow  in 
the  last  section. 


2  High  order  Roe  scheme  for  nonconservative  hyper¬ 
bolic  systems 

In  this  section,  we  follow  the  procedure  in  [7]  to  define  the  high  order  Roe  scheme  for 
solving  nonconservative  hyperbolic  systems.  For  the  one-dimensional  case,  a  nonconservative 
hyperbolic  system  reads 

Wt  +  A(W)WX  =  0,  xen  CM,  t>  0  (2.1) 

where  W  =  W{x,  t)  is  a  N-component  state  vector  and  A(W)  is  a  N  x  N  matrix.  The  system 
is  supposed  to  be  hyperbolic,  i.e.  A(W)  has  N  real  eigenvalues  and  a  full  set  of  N  linearly 
independent  eigenvectors. 

We  use  a  uniform  grid 

a  =  xi  <  X3  <  ■  ■  ■  <  xN  i  <  xN  ,  i  =  b. 

2  2  2  iVa:+2 

The  cells,  cell  centers,  and  the  uniform  cell  size  are  denoted  by 

h  =  [xi_ i,xi+i],  Xi  =  ^(Xi_ i+xi+i),  Ax  =  xi+i  i  =  l,2,...,Nx. 

In  the  case  of  systems  of  conservation  laws,  that  is  when  A(W)  =  dF/dW,  which  is  the 
Jacobian  of  a  flux  function  F(W),  (2.1)  reduces  to  a  classical  conservation  law 

Wt  +  F(W)x  =  0  (2.2) 

A  conservative  hnite  volume  semi-discretization  of  the  system  (2.2)  is 

Wl(t)  =  ^(Gj-i/2  —  Gi+ 1/2),  (2.3) 

with  the  numerical  flux 

Gi+i/ 2  =  G(W~1/2(t),  W+1/2(t)),  (2.4) 
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where  Wt(t)  is  used  to  approximate  the  cell  averaged  value  Wi(t),  which  is  dehned  as 


_  i  rxi+ 1/2 

Wi(t)  =  —  /  W(x,t)dx, 

Ax  4-1/2 

and  H/^1/2(t)  are  the  reconstructed  states  associated  to  the  cell  average  sequence 
The  semi-discrete  high  order  Roe  scheme  for  (2.3)  can  be  written  as 

W7M  =  -  w~1/2(t)) 

+A+1  ~  w/;i/2 (<))  -  ^(^-1/2(0)  +  m,;1/2W))  (2.5) 

which  is  equivalent  to  the  conservative  scheme  (2.3)  with  the  numerical  flux  (2.4)  dehned  to 
be 

Gi+1/2  =  F(W-+in(t))  +  A7+in(W++1/2(t)  -  WMn(t)) 

=  -nW*+ll2(t))  +  At+in(W*ll2(t)-WMI2(t))  (2.6) 

=  +  F(Ww/2(t ))  - \Awn\(W+ll2(t)  -  ww/2(t))). 

Here 

|4+1/2|  =  4+1/2  _  4+1/2  (2-7) 

and  the  Roe  property 

A(W+  -  W~)  =  F(W+)  -  F(W~)  (2.8) 

has  been  used.  The  intermediate  matrices  are  dehned  by 

A+i/2  =  +wy1/2(f),  wq1/2(f))  (2.9) 

and 

4+1/2  —  ^*+1/24+l/24+l/2’  4+1/2  —  di(lQ  ((4)j_|_i/25  )  4a4+1/2)  (2-10) 

where  R4+1/2  is  a  N  x  N  matrix  with  each  column  as  a  right  eigenvector  of  Ai+ 1/2,  and  Aj+i/2 
is  the  diagonal  matrix  whose  diagonal  entries  are  the  corresponding  eigenvalues  of  Ai+i/2- 
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We  introduce  P/(x)  as  any  smooth  function  defined  in  the  cell  such  that 

lim  P‘(x)  =  W+1/2(i),  lim  P‘{x)  =  W,;1/2(«).  (2.11) 

X^Xi-l/2  x^xi+ 1/2 

Then,  (2.5)  can  be  rewritten  as 

+AM/2(W+1/2(t)  -  W-+1/2(t))  +  /  F(P‘(x))dx).  (2.12) 

Xxi-l/2  UX 

Now  the  numerical  high  order  Roe  scheme  for  solving  the  nonconservative  system  (2.1) 
can  be  easily  generalized  from  (2.12) 

WRt)  =  -^(At1/2(W+1/2(t)-W-1/2(t)) 

rxi+i/2  n 

+A_+./2(M/ii/2(‘)  -  Wrn^W)  +  /  A(P‘(x))-P‘(x)dx)  (2.13) 

Jxi_x/2  aX 

with  the  function  P*(x)  satisfying  (2.11).  In  order  to  obtain  entropy-satisfying  solutions,  the 
Harten-Hyman  entropy  fix  technique  [17,  44]  can  be  applied  to  this  scheme. 

3  Two-medium  flow  model  for  nonconservative  Euler 
equations 

In  this  section,  we  describe  the  two-medium  inviscid  compressible  flow  model  for  the  non¬ 
conservative  or  primitive  Euler  equations.  As  pointed  out  in  [21],  the  choice  of  the  primitive 
set  of  variables,  that  include  density,  velocity  and  pressure,  provides  a  model  better  suited 
for  computations  of  propagating  material  fronts  and  results  in  clean  and  monotonic  solution 
profiles,  so  we  consider  the  one  dimensional  primitive  Euler  equations 

Wt  +  A(W)WX  =  0  (3.1) 

with 

(up  0  \ 

W  =  (p,u,p)T,  A(W)=  0  u  p-1  .  (3.2) 

yo  7 p  u  ) 
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Here  p  is  the  density,  u  is  the  velocity,  p  is  the  pressure,  7  is  the  ratio  of  specific  heats.  The 
total  energy  is  given  by  E  =  pe  +  \pu 2,  where  e  is  the  specific  internal  energy  per  unit  mass. 
The  7-law  equation  of  state  (EOS)  used  for  ideal  gases  is  given  as 

pe=p/(  7-1)  (3-3) 

and  Taft  EOS  used  for  the  water  medium  [9,  13,  25,  33]  is  expressed  as 

pe=(p  +  NwB)/{Nw-  1)  (3.4) 

where  B  =  B  —  A,  Nw  =  7.15,  A  =  1.0E5  Pa  and  B  =  3.31E8  Pa. 

4  WENO  reconstruction  with  subcell  resolution 

In  this  section,  we  will  describe  how  to  use  the  WENO  reconstruction  with  subcell  resolu¬ 
tion  to  reconstruct  from  the  cell  averages  The  WENO  reconstruction  is 

described  in  detail  in  [20,  37].  We  follow  the  procedure  in  [39]  to  describe  how  to  apply  the 
subccll  resolution  technique  of  Harten  [15]  to  the  scheme  (2.13)  with  the  third  order  TVD 
Runge-Kutta  time  discretization  [38]. 

We  first  consider  the  ID,  scalar,  linear  version  ut  +  f{u)x  =  0,  with  f(u )  =  au  and  a  >  0, 
to  describe  the  WENO  reconstruction  with  the  subcell  resolution  technique.  The  extension 
to  the  nonlinear  and  system  cases  will  follow.  We  would  like  to  reconstruct  u~+1,2  and  u 4w2 
in  each  cell  from  the  sequence  of  cell  averages  {77,},  with  the  following  algorithm. 

WENO  Reconstruction  Algorithm: 

Given  the  cell  averages  {77^}  of  a  function  u(x ): 

1  rxi+ 1/2 

uj  =  —  u(£)d£,  j  =  1,2,...,  Nx,  (4.1) 

^  ^-1/2 

based  on  the  big  stencil  Si  =  {Ii-r,  ■  ■  ■  ,  •  •  •  ,  Ii+r},  a  fc-th  (k  =  2r  +  1)  order  accurate 

approximation  to  the  boundary  values  u~+lj2  and  2  in  the  cell  is  reconstructed  as 

r  r 

uI+i/2  =  J2^j(xi+i/2)pj(xi+1/2),  U+_  1/2  =  'Y^Uj(xi-x/2)Pj{Xi-l/-l)  (4.2) 

7=0  7=0 

8 


where  each  p3  (x)  is  a  reconstruction  polynomial  that  uses  the  cell  averages  in  the  small 
stencil  Sf  =  {h-j,  •  •  •  ,Ii-j+r}  C  S).  The  nonlinear  weights  {ujj(x)}r3=0  are  calculated  from 
the  polynomials  {Pj(x)}r3=0  and  the  linear  weights  {dj(x)}rj=0  at  each  fixed  point  x,  and  they 
satisfy 

r 

Uj(x)>  o,  ^cjj(x)  =  1.  (4.3) 

3=0 

The  nonlinear  weight  uj3(x)  is  close  to  zero  when  a  discontinuity  is  located  in  the  stencil  S3, 
so  as  to  avoid  involving  much  information  from  any  stencil  Sf  which  contains  discontinuities. 

Remark:  For  the  cell  boundaries  (4.2),  the  linear  weights  d3{xl+\/2)  and  dj{xi- 1/2)  are 
positive.  However,  at  certain  internal  points  x  e  R  (reconstruction  at  those  points  are 
needed  in  Section  5),  the  linear  weights  d3  (x)  may  be  negative.  The  linear  weights  may  also 
become  negative  if  the  stencil  Si  is  changed  to  Si+1  or  St_\ ,  while  still  reconstructing  values 
in  the  cell  R  (e.g.  in  the  following  subcell  resolution  algorithm).  In  these  cases,  the  technique 
to  treat  negative  weights  in  [36]  needs  to  be  applied. 


Subcell  Resolution  Algorithm: 

First,  at  the  beginning  of  every  Runge-Kutta  cycle: 

(I)  Define  the  “critical  intervals”  (intervals  containing  discontinuities)  /,;  =  (xj_i/2,  £*+1/2) 
by  cq  >  ai+1,  Gi  >  CTj-i,  where  at  =  \m(/\+uh  A_ttj)|,  A +Ui  =  ui+1  -  ui:  A =  u{  -u^ x 
and  m  is  the  minmod  function  which  is  defined  to  be 


,  •  •  •  ,  On) 


s  minx<i<n  1^1^ 

0, 


if  s  =  sign  (a  1)  =  •  •  •  =  sign(an) 
otherwise  . 


(4.4) 


(II)  For  any  “critical  interval”  II.  let  6*j  =  (ui+ 1  —  Ui)/(ui+ 1  —  1),  and  use  Xi_i/2  +  djAx 

as  an  approximation  to  the  discontinuity  location  inside  the  cell  It. 

Then,  in  each  Runge-Kutta  time  cycle,  we  perform: 

(III)  Let  the  cell  I,,  boundary  values  u~+1/2  and  be  dehned  as  usual,  using  the 

standard  WENO  reconstruction  algorithm  (4.2),  unless  It  or  (for  the  second  and  third  Runge- 
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Kutta  stages)  It-i  is  a  “critical  interval”.  If  Jj  is  a  “critical  interval”,  we  define 


Ui+ 1/2 


(4.5) 


—  .old 
Ui- 1/2 


(4.6) 


where  is  the  standard  WENO  reconstruction  with  the  stencil  Si- i,  m,+1/2  is  also  the 

standard  WENO  reconstruction  with  the  stencil  Si- 1,  but  evaluated  at  xi+i/2  of  the  cell  /*, 


( Ft} 

and  u)+y2  is  the  WENO  reconstruction  with  the  stencil  Si+±  and  evaluated  at  Xi+ \/2.  Notice 
that  here  for  W-+1/2  the  technique  for  treating  negative  weights  needs  to  be  used.  For  the 
second  and  third  Runge-Kutta  stages,  we  choose  the  stencil  S)+ 2  for  u\+]/2  h  0  <  1  and  the 
negative  weight  treating  technique  needs  to  be  used  here  as  well;  and  when  /*_!  is  a  “critical 


interval”  and  <  1,  we  choose  the  stencil  Si+ 1  for  ui+1^2. 

Remark: 

(a)  The  case  for  a  <  0  is  easily  obtained  by  symmetry. 

(b)  If  two  adjacent  cells  are  both  “critical  intervals”,  then  we  remove  the  one  with  a 
smaller  cq  from  the  list  of  critical  intervals. 

(c)  For  nonlinear  systems,  the  subcell  resolution  algorithm  is  simply  applied  in  each  local 

characteristic  held.  The  detailed  procedure  can  be  found  in  [39].  The  only  difference  in  our 
current  situation  is  that  we  just  use  A{Wi)  as  the  Roe  average  matrix  to  do  the  characteristic 
decomposition  when  defining  if  is  a  “critical  interval” .  However  if  the  material  interface  is 
located  in  the  cell  Jj,  we  always  set  it  to  be  a  “critical  interval”,  and  we  use  A  to 

do  the  characteristic  decomposition  for  the  right  cell  boundary  at  £*+1/2,  and  A  j 

for  the  left  cell  boundary  at  \/2.  The  Roe  average  to  define  the  intermediate  matrix  in 


(d)  It  is  pointed  out  in  [15,  39]  that  the  subcell  resolution  technique  should  be  applied 
only  to  sharpen  contact  discontinuities.  Special  caution  is  needed  when  one  tries  to  sharpen 
a  (nonlinear)  shock,  to  avoid  obtaining  a  nonphysical,  entropy  condition  violating  solution. 
In  our  approach,  the  subccll  resolution  is  used  to  get  sharp  left  and  right  states  at  the 
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discontinuities.  At  a  nonlinear  shock,  this  is  also  needed  so  as  to  get  a  more  accurate  shock 
speed  which  heavily  depends  on  the  left  and  right  states.  In  the  computation  for  Euler 
equations  of  compressible  gas  dynamics,  we  apply  the  subcell  resolution  algorithm  in  both  the 
linearly  degenerate  held  and  genuinely  nonlinear  fields  [44],  but  for  the  genuinely  nonlinear 
fields,  with  eigenvalues  \L  and  XR  corresponding  to  the  left  and  right  states  respectively, 
we  only  apply  the  subcell  resolution  for  the  case  AL  >  \R  to  avoid  sharpening  a  rarefaction 
wave. 


5  Choice  of  the  path  and  evaluation  of  the  path  inte¬ 
gral 


In  smooth  regions,  all  simple  wave  models  for  conservative  Euler  equations  and  primitive 
Euler  equations  are  equivalent.  However,  near  discontinuities,  they  are  not  equivalent  [21]. 
According  to  this,  our  choice  of  the  path  is  divided  into  three  parts:  the  smooth  case, 
discontinuities  in  a  single  medium,  and  discontinuities  at  the  material  interface. 

5.1  The  smooth  case 


In  the  smooth  case,  the  integral  term  in  (2.13) 


fxi+l/2  (I 

/  A(P*(x))-Pt(x)dx 

'xi- 1/2  aX 


(5.1) 


can  be  computed  via  a  high  order  accurate  Gauss-Lobatto  quadrature  rule.  Given  the 
positions  {sj}  and  associated  weights  {u>j}  for  a  G-point  quadrature  in  the  interval  [— |], 
we  can  replace  the  analytical  path  integral  (5.1)  by 

G 

,  Pt(  T.\rl.T.  = 

dx 


xi+ 1/2  r]  -  (j 

A(P‘(x))—P‘(x)dx  =  '£,“jMp!(sj))—P‘(sj). 

i- 1/2  j  =  1 


(5.2) 


In  our  numerical  experiments,  we  use  the  four-point  Gauss-Lobatto  quadrature  rule: 


Sl’4_T2’ 


S2,3  =  =F' 


,V5 
10  ’ 


^1,4  — 


<^2,3  — 


(5.3) 
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For  Pt(s1A),  from  (2.11)  we  already  have 

P?M  =  W+1/2,  P‘(s  4)  =  W~+1/1  (5.4) 

following  the  same  procedure  for  obtaining  P/(s i)4),  we  can  also  obtain  P/(s2i 3),  where  P/(s2) 
is  obtained  in  the  same  way  as  P\ (s  1 )  corresponding  to  the  local  characteristic  field  at  Xj_4/2, 
and  P*{s3)  as  P/(s4)  corresponding  to  the  local  characteristic  field  at  xi+i/2.  Note  here  for 
the  smooth  case,  we  do  not  need  the  subccll  resolution  in  the  WENO  reconstruction.  Since 
we  have  {P/(sj)}^=  4,  ^P/(x)  can  be  approximated  by  Q{x),  the  derivative  of  the  Lagrangian 
interpolation  polynomial  based  on  {P/(s\,)}d=1.  Then  -^P^(sj)  can  be  replaced  by  Q(sj)  in 
(5.2). 

5.2  Discontinuities  in  a  single  medium 

At  the  discontinuities,  the  cell  is  defined  as  a  “critical  interval” ,  and  we  have  obtained  the  left 
and  right  states  W+_ , and  W~+1j2  from  the  WENO  reconstruction  with  subcell  resolution 
in  Section  4.  Denote  the  left  and  right  states  W+_l/2  and  ,2  to  be  Wl  and  Wr,  we  can 
use  the  exact  Riemann  solver  [40,  44]  to  obtain  the  exact  Riemann  solution  between  the 
two  states.  The  exact  Riemann  solution  for  the  compressible  Euler  equations  contains  four 
constant  states  connected  by  a  rarefaction  wave  or  shock  wave,  a  contact  discontinuity,  and 
another  rarefaction  wave  or  shock  wave.  The  four  constant  states  can  be  denoted  as  Wl, 
W*l,  W*r  and  Wr.  In  this  case  of  the  “critical  interval”,  we  use  the  exact  Riemann  solution 
to  get  the  integral  path  for  the  integral  term  (5.1).  It  can  be  computed  as  the  following 
seven  parts: 

(a)  The  four  constant-state  parts,  since  j-P/(x)  =  0,  the  integral  in  these  four  parts  are 
all  zero. 

(b)  If  Wl  and  W*l  are  connected  by  a  rarefaction  wave,  we  similarly  use  the  4-point 
Gauss-Lobatto  quadrature  rule  (5.3),  and  P/(x)  in  this  part  is  just  the  rarefaction  wave  line. 
Otherwise,  if  Wl  and  W*l  are  connected  by  a  shock  wave,  then  the  integral  path  for  this  part 
needs  to  satisfy  the  Rankine-Hugoniot  jump  condition  of  the  conservative  Euler  equations, 


12 


and  the  integral  result  is  o{W*l  —  Wl),  with  a  being  the  shock  speed  related  to  the  two 
states  Wl  and  W*l-  Similar  results  can  be  obtained  for  the  connection  between  W*r  and 
WR. 

(c)  For  the  contact  discontinuity  part  between  W*l  and  W*r,  the  integral  path  also 
needs  to  satisfy  the  Rankine-Hugoniot  jump  condition.  Since  W*l  and  have  the  same 
pressure  p  and  velocity  u  and  different  densities  p*L  and  p*R,  the  integral  result  is  simply 
-  P*l),  0,0)t. 

Remark:  In  (b)  and  (c)  above,  for  the  shock  wave  and  contact  discontinuity,  we  do  not 
need  to  know  the  exact  integral  paths  for  satisfying  the  Rankine-Hugoniot  jump  condition, 
in  order  to  get  the  integral  results  along  those  paths.  We  only  need  to  make  sure  that  these 
paths  exist,  which  can  be  easily  verified. 

5.3  Discontinuities  at  the  material  interface 

We  always  set  the  cell  at  the  material  interface  as  a  “critical  interval”.  The  integral  term 
(5.1)  can  be  computed  similar  to  the  second  case  of  discontinuities  in  a  single  medium,  as 
the  exact  Riemann  solution  at  the  material  interface  is  almost  the  same  as  that  for  the  single 
component  compressible  Euler  equations,  also  containing  four  constant  states  connected  by 
a  rarefaction  wave  or  a  shock  wave,  a  contact  discontinuity,  and  another  rarefaction  wave 
or  shock  wave  [25,  27,  24],  Apart  from  W*l  and  W*r  which  are  connected  by  a  contact 
discontinuity,  the  left  part  is  for  medium  one  with  the  ratio  of  specific  heats  71  and  the  right 
part  is  for  medium  two  with  the  ratio  of  specific  heats  72. 

6  Tracking  the  moving  medium  interface 

In  this  section,  we  describe  how  to  use  the  level  set  equation  [41,  4,  31]  to  track  the  moving 
fluid  interface.  The  level  set  equation  for  the  one- dimensional  case  is 

<f>t  +  u<f>x  =  0.  (6.1) 
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The  interface  is  tracked  as  the  zero  level  set  of  0,  with  the  initialized  0(x)  to  be  the  signed 
normal  distance  to  the  material  front.  We  use  a  fifth  order  finite  difference  WENO  method 
[20,  37]  with  the  third  order  TVD  Runge-Kutta  time  discretization  [38]  to  solve  the  level 
set  equation  (6.1).  This  equation  is  solved  concurrently  with  the  nonconservative  Euler 
equations  (3.1),  using  the  velocity  u  coming  from  the  Euler  equations.  The  solution  0  =  0 o 
from  solving  the  level  set  equation  (6.1)  has  the  zero  level  set  as  the  material  interface,  but 
it  needs  not  be  the  distance  function  for  t  >  0.  A  serious  distorted  level  set  function  0  =  0O 
may  lead  to  significant  errors  for  t  >  0.  For  this  reason,  0 o(x)  is  reinitialized  to  be  a  signed 
normal  distance  function  to  the  interface  by  solving  the  following  eikonal  equation  to  steady 
state 

0r  =  5(00) (1  -  |0x|)  (6.2) 

through  iterating  the  pseudo-time  r,  where  =  -j==t==  is  the  approximate  sign  func¬ 
tion.  We  use  the  Godunov  Hamiltonian  and  the  fifth  order  finite  difference  WENO  discretiza¬ 
tion  in  [19,  50]  to  solve  (6.2).  The  stopping  criterion  for  this  iteration  is  e\  <  Ar(Ax)2,  where 
ei  is  the  L 1  difference  of  0  between  two  consecutive  iteration  steps,  and  we  take  At  =  Ax/10 
in  the  experiments  [41]. 

7  Algorithm  summary 

Our  basic  semi-discrete  scheme  is  (2.13),  which  can  be  written  as 

Wt  =  L(W). 

It  is  discretized  in  time  by  the  third  order  TVD  Runge-Kutta  method  [38]: 


=  Wn  +  AtL(Wn), 

w{2) 

=  - Wn  +  ifK(1)  +  -A  tL(Ww), 

4  4  4 

(7.1) 

wn+1 

=  +  ^A  t.L(W{2)). 

o  o  o 

Notice  that  this  Runge-Kutta  method  is  simply  a  convex  combination  of  three  Euler  forwards. 
We  can  now  summarize  our  algorithm  to  advance  one  Euler  forward  time  step  (still  denoted 
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as  from  tn  to  tn+1)  in  the  following  steps. 


Step  1.  Compute  the  new  time  step  size  based  on  the  CFL  condition: 

At  =  CFL  Ax/  max  (|<|  +  c?)  (7.2) 

1  <j<Nx  J  J 

where  c”  =  sj^fp/ / p™  is  the  sound  speed,  and  p™,  u/\p/  are  the  density,  velocity  and  pressure 
at  time  level  tn,  respectively.  This  step  needs  to  be  done  only  at  the  beginning  of  the  whole 
Runge-Kutta  cycle. 

Step  2.  Taking  Wn  and  o'1  as  the  initial  condition,  solve 


Wt  =  L(W), 


(7.3) 


<h  =  p(4>). 


(7.4) 


for  one  time  step  using  the  Runge-Kutta  time  discretization  (7.1).  Here  (7.4)  is  written  from 
(6.1)  with  the  fifth  order  finite  difference  WENO  spatial  discretization.  At  each  Runge- 
Kutta  time  cycle,  we  first  reconstruct  and  WVW2  as  described  in  Section  4,  then, 

based  on  W/r_x,2  and  W~+1j2,  we  compute  the  integral  term  (5.1)  as  described  in  Section  5. 
We  can  then  formulate  the  right  side  of  (7.3).  P(«p)  can  be  formulated  simultaneously,  given 
the  velocity  u  =  un  from  the  Euler  equations.  Denote  the  updated  W  by  kFn+1,  and  the 
updated  0  by  0n+1/2. 

Step  3.  Reinitialize  0n+1/2  by  solving  (6.2)  with  0o  =  0n+1//2,  and  denote  this  solution  by 
(j)n+1. 

Step  4.  Define  the  new  interface  position  from  the  zero  level  set  of  (pn+1 .  Now  we  have 
advanced  one  Euler  forward  time  step. 


8  Numerical  experiments 

In  the  following,  we  show  several  one- dimensional  numerical  examples  to  demonstrate  that 
our  approach  can  improve  the  performance  of  the  high  order  Roe  scheme  for  the  noncon¬ 
servative  Euler  equations.  With  fifth  order  WENO  reconstruction  in  Section  4  (r  =  2)  and 
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four-point  Gauss-Lobatto  quadrature  rule  in  Section  5.1,  the  Roe  scheme  can  achieve  fifth 
order  accuracy  for  smooth  solutions,  which  will  be  tested  in  Example  1.  Except  for  Example 
1,  Example  3  and  Example  5,  the  computational  domain  for  all  other  examples  are  taken 
as  [0,1],  and  the  initial  material  interface  for  the  two-medium  flow  problems  is  located  at 
x  =  0.5.  The  CFL  number  is  taken  as  0.1.  The  computational  domain  is  divided  with 
Nx  =  100  uniform  grids.  Units  for  density,  velocity,  pressure,  length  and  time  are  kg/m 3, 
m/s,  Pa,  m  and  s,  respectively. 

8.1  Example  1 

In  this  example,  we  first  test  the  accuracy  for  a  smooth  solution  to  a  single  component 
nonconservative  Euler  equations  with  initial  conditions 

Po(x )  =  1  +  0.8  sin(a:),  u0(x )  =  1,  Po{x)  =  1  (8.1) 

on  a  domain  [0,  2ir\  with  periodic  boundary  conditions.  The  exact  solution  is 

p{x)  =  1  +  0.8  sin(x  —  t ),  u{x )  =  1,  p[x )  =  1. 


Since  the  solution  to  this  problem  is  smooth,  we  do  not  need  to  apply  the  subccll  resolution. 
With  fifth  order  WENO  reconstruction  in  Section  4  (r  =  2)  and  four-point  Gauss-Lobatto 
quadrature  rule  for  the  integral  term  in  the  smooth  case  in  Section  5.1,  the  fifth  order 
accuracy  can  be  achieved  for  this  path-conservative  scheme  applied  to  the  single  component 
nonconservative  Euler  equations,  as  listed  in  Table  8.1. 


Table  8.1:  Accuracy  test  for  Example  1  with  initial  data  (8.1),  t=l. 


Nx 

Ll  error 

order 

L°°  error 

order 

20 

2.02E-04 

- 

3.74E-04 

- 

40 

5.96E-06 

5.08 

1.28E-05 

4.87 

80 

1.81E-07 

5.04 

3.93E-07 

5.03 

160 

5.58E-09 

5.02 

1.19E-08 

5.05 

320 

1.76E-10 

4.99 

3.53E-10 

5.07 
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8.2  Example  2 


In  this  example,  we  now  show  a  single  medium  flow  problem,  which  is  the  standard  Lax 
shock  tube  problem,  to  demonstrate  that  we  obtain  the  correct  entropy  solution  by  our 
non-conservative  scheme.  Here  7  =  1.4  is  used,  with  the  initial  condition: 

(0.445,0.698,3.528)  for  x  <  0.5; 


(P,  P )  = 


(8.2) 


(0.5,0,0.571)  for  x  >  0.5. 

The  computed  density  p,  velocity  u  and  pressure  p  are  plotted  at  t  =  0.17  against  the  exact 
solution  in  Fig.  8.1.  We  can  see  that  the  rarefaction  wave,  contact  discontinuity  and  shock 
wave  are  all  captured  well. 

8.3  Example  3 

This  example  is  a  right  moving  shock  for  a  single  medium  flow.  Here  7  =  5/3  is  used,  with 
the  initial  condition: 

;fi,9y^,10)  for  x  <  0.5; 

(1,0,1)  for  x  >  0.5. 

We  compute  this  problem  on  a  domain  [0,1].  A  similar  example  (in  Lagrangian  form)  is 
used  in  [3]  to  demonstrate  potential  problems  of  nonconservative  path-based  Roe-type  and 
Lax-Friedrichs  (LxF)-type  schemes.  The  computed  density  p,  velocity  u  and  pressure  p  are 
plotted  at  t  —  0.1  against  the  exact  solution  in  Fig.  8.2.  We  can  see  that  our  algorithm  can 
capture  the  correct  shock  location  in  this  case,  and  the  minor  glitches  on  the  left  constant 
state  would  not  grow  with  a  much  refined  mesh  (figures  on  the  right). 

We  also  use  this  simple  example  to  illustrate  that  in  our  approach,  the  subccll  resolution 
and  the  exact  Riemann  solution  are  both  necessary  in  order  to  catch  the  right  path  in  this 
non-conservative  scheme.  In  Fig.  8.3,  we  list  the  results  for  two  cases:  one  is  obtained  when 
we  use  the  exact  Riemann  solution  but  without  subcell  resolution,  the  other  is  obtained 
when  we  use  subcell  resolution  but  with  the  line  path  [43,  3]  instead  of  the  exact  Riemann 
solution,  on  a  very  refined  mesh  Nx  =  1000.  We  can  see  that  neither  of  them  can  catch  the 
correct  shock  solution. 
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1.4 


Figure  8.1:  Density,  velocity  and  pressure  for  Example  2.  t  =  0.17.  Solid  line:  the  exact 
solution.  Symbol:  the  numerical  solution. 
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Figure  8.2:  Density,  velocity  and  pressure  for  Example  3.  t  =  0.1.  Solid  line:  the  exact 
solution.  Symbol:  the  numerical  solution.  Left:  mesh  Nx  =  100;  Right:  mesh  Nx  =  1000. 
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3.5 


3.5 


0.2  0.4  0.6  0.8 


0.2  0.4  0.6  0.8 


Figure  8.3:  Density,  velocity  and  pressure  for  Example  3.  t  =  0.1.  Nx  =  1000.  Solid  line: 
the  exact  solution.  Symbol:  the  numerical  solution.  Left:  exact  Riemann  solution  without 
subcell  resolution;  Right:  subcell  resolution  and  with  the  line  path  instead  of  the  exact 
Riemann  solution. 
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8.4  Example  4 


This  is  an  air-helium  shock  tube  problem  taken  from  [13,  33,  28],  with  the  initial  condition: 

(8.4) 


s  _  i  (1,  0, 1  x  105, 1.4) 
ip,u,P,l)  ^  (0.125,0,1  x  104, 1.2) 


for  x  <  0.5; 
for  x  >  0.5. 


The  computed  density  p,  velocity  u  and  pressure  p  are  plotted  at  t  —  0.0007  against  the 
exact  solution  in  Fig.  8.4.  We  can  see  that  the  rarefaction  wave,  the  contact  discontinuity 
and  the  shock  wave  are  all  captured  well  for  this  two-medium  flow. 

8.5  Example  5 

This  example  is  used  to  demonstrate  the  advantage  of  high  order  methods  as  in  [47]  for 
two  medium  flows.  It  contains  both  shocks  and  fine  structures  in  smooth  regions,  which  is 
a  simple  model  for  shock-turbulence  interactions.  The  initial  discontinuity  is  right  on  the 
material  interface  and  located  at  x  =  —4,  the  left  and  right  states  for  the  initial  discontinuity 
are 


(P,u,p,  7)  = 


(3.857143,  2.629369, 10.333333, 1.4) 


for  x  <  —4, 
for  x  >  —4. 


(8.5) 


(1  +  0.2  sin(5x),  0, 1, 5/3) 

We  compute  the  example  up  to  time  f  =  1.8  on  a  domain  [—5,  5].  For  comparison,  we  also 
show  the  results  obtained  with  a  second  order  method,  by  replacing  the  fifth  order  WENO 
reconstruction  with  subccll  resolution  in  Section  4  with  a  second  order  ENO  reconstruction 
[15,  39],  and  replacing  the  four-point  Gauss-Lobatto  quadrature  rule  for  the  smooth  integral 
term  in  Section  5.1  with  a  trapezoid  rule.  We  plot  the  densities  obtained  by  the  second 
order  and  the  fifth  order  methods  with  Nx  =  200  in  Fig.  8.5.  The  solid  line  is  the  solution 
obtained  by  the  fifth  order  method  with  Nx  =  2000  points,  which  can  be  considered  as 
a  converged  reference  solution.  We  can  find  that  the  two  methods  can  both  capture  the 
correct  solution,  however  the  result  of  the  fifth  order  method  is  in  better  agreement  with 
the  converged  reference  solution  than  the  second  order  method,  especially  in  the  region  with 
fine  structures,  on  this  relatively  coarse  mesh.  This  is  is  similar  to  the  results  in  [47]. 
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1.2 


Figure  8.4:  Density,  velocity  and  pressure  for  Example  4.  t  =  0.0007.  Solid  line:  the  exact 
solution.  Symbol:  the  numerical  solution. 
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Figure  8.5:  Density,  velocity  and  pressure  for  Example  5.  t  =  1.  Solid  line:  the  converged 
reference  solution  for  the  fifth  order  method  with  Nx  =  2000.  Symbol:  the  numerical  solution 
with  Nx  =  200.  Left:  second  order  method;  Right:  fifth  order  method. 

8.6  Example  6 

This  is  a  problem  of  a  shock  wave  refracting  at  an  air-helium  interface  with  a  reflected  weak 

rarefaction  wave  taken  from  [13,  33,  28],  with  the  initial  condition: 

{(1.3333,  0.3535\/l05,  1.5  x  105, 1.4)  for  x  <  0.05, 

(1,0,1  x  105, 1.4)  for  0.05  <x<  0.5,  (8.6) 

(0.1379,0,1  x  105,  5/3)  for  x  >  0.5. 

The  computed  density  p,  velocity  u  and  pressure  p  are  plotted  at  t  —  0.0012  against  the 

exact  solution  in  Fig.  8.6.  The  strength  of  the  shock  for  this  example  is  Pl/pr  =  1.5, 

the  computed  results  compare  well  with  the  exact  solutions,  without  oscillation  around  the 

material  interface  for  the  density. 


8.7  Example  7 

This  example  is  the  same  as  Example  6,  also  taken  from  [13,  33,  28]  only  by  increasing  the 

strength  of  the  right  shock  wave  to  Pl/pr  =  15,  with  the  initial  condition  as: 

{(4.3333,  3.2817-^105, 1.5  x  106, 1.4)  for  x  <  0.05, 

(1,0,1  x  105, 1.4)  for  0.05  <  x  <  0.5,  (8.7) 

(0.1379,0,1  x  105,  5/3)  for  x  >  0.5. 

The  computed  density  p,  velocity  u  and  pressure  p  are  plotted  at  t  —  0.0005  against  the 

exact  solution  in  Fig.  8.7.  The  computed  results  still  compare  reasonably  well  to  the  exact 
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Figure  8.6:  Density,  velocity  and  pressure  for  Example  6.  t  =  0.0012.  Solid  line:  the  exact 
solution.  Symbol:  the  numerical  solution. 
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solutions,  and  there  is  still  no  oscillation  around  the  material  interface  for  the  density  with 
this  stronger  strength  of  the  shock.  The  small  glitches  on  the  left  constant  state  can  also 
be  observed  for  the  methods  in  [13,  33,  28],  which  were  explained  in  [13]  to  be  due  to  the 
(mis)capturing  of  the  perfect  shock  initial  data  by  a  shock  capturing  scheme.  In  this  example, 
it  is  more  pronounced  since  the  shock  wave  is  much  stronger  compared  to  Example  6. 

8.8  Example  8 

This  is  a  problem  of  a  shock  wave  refracting  at  an  air-helium  interface,  but  with  a  reflected 
shock  wave  [13,  33],  with  the  initial  condition  as: 


( (1.3333,  O^SS^/lO^,  1.5  x  105, 1.4) 

for 

x  <  0.05, 

(p,  U,  p,  7)  =  < 

(1,0,1  xlO5, 1.4) 

for 

0.05  <  x  <  0.5, 

(8.8) 

[(3.1538,0,1  x  105, 1.249) 

for 

x  >  0.5. 

The  computed  density  p,  velocity  u  and  pressure  p  are  plotted  at  t  —  0.0017  against  the 
exact  solution  in  Fig.  8.8.  The  computed  results  show  that  there  are  slight  oscillations  near 
the  left  shock  region,  similar  results  can  also  be  observed  in  [13]  using  the  standard  scheme 
from  [30].  If  we  increase  the  strength  of  the  shock  to  Pl/pr  =  15,  the  code  would  blow 
up.  It  appears  that  our  subcell  resolution  procedure  is  still  not  accurate  enough  to  capture 
the  strong  generated  shock  during  the  instant  when  the  original  shock  wave  impacts  on  the 
material  interface. 


8.9  Example  9 


This  is  a  gas-water  shock  tube  problem  with  very  high  pressure  in  the  gaseous  medium  taken 
from  [33,  28].  The  initial  condition  is  given  as 


(P,u,p,  7)  = 


for  x  <  0.5; 
for  x  >  0.5. 


(8.9) 


I  (1270,0,8  x  108, 1.4) 

\  (1000,  0, 1  x  105,  7.15) 

The  computed  density  p,  velocity  u  and  pressure  p  are  plotted  at  t  —  0.00016  against  the 

exact  solution  in  Fig.  8.9.  Even  though  the  initial  pressure  in  the  gas  is  extremely  high  and 

a  very  strong  shock  is  generated  in  the  water,  our  computed  results  still  compare  well  to  the 

exact  solutions. 
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Figure  8.7:  Density,  velocity  and  pressure  for  Example  7.  t  =  0.0005.  Solid  line:  the  exact 
solution.  Symbol:  the  numerical  solution. 
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Figure  8.8:  Density,  velocity  and  pressure  for  Example  8.  t  =  0.0017.  Solid  line:  the  exact 
solution.  Symbol:  the  numerical  solution. 
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Figure  8.9:  Density,  velocity  and  pressure  for  Example  9.  t  =  0.00016.  Solid  line:  the  exact 
solution.  Symbol:  the  numerical  solution. 


8.10  Example  10 


This  example  increases  the  energy  of  the  explosive  gaseous  medium  in  Example  9,  which 
is  taken  from  [13,  33,  28],  originally  studied  in  [49]  for  the  underwater  explosions,  with  the 
initial  condition  given  as 


for  x  <  0.5; 
for  x  >  0.5. 


(8.10) 


The  computed  density  p,  velocity  u  and  pressure  p  are  plotted  at  t  —  0.0001  against  the 
exact  solution  in  Fig.  8.10.  We  begin  to  see  some  discrepancies  between  the  exact  solution 
and  the  numerical  solution,  however  the  errors  are  still  reasonably  small  considering  that  we 
are  using  a  nonconservative  scheme  on  this  example  with  very  strong  discontinuities.  The 
methods  in  [13,  33,  28]  can  solve  this  problem  well. 

9  Concluding  remarks 

In  this  paper,  we  have  investigated  using  the  WENO  reconstruction  with  a  subcell  resolution 
technique  when  applied  to  high  order  Roe-type  schemes  for  nonconservative  Euler  equations. 
The  subcell  resolution  to  sharpen  discontinuities,  in  order  to  remove  or  significantly  reduce 
transition  points  at  discontinuities,  has  been  shown  to  significantly  improve  the  capturing 
of  the  correct  discontinuity  waves  by  defining  the  correct  path  in  the  path  integral  at  dis¬ 
continuities.  While  this  paper  identifies  the  smearing  of  discontinuities  by  shock  capturing 
schemes  with  transition  points  at  the  discontinuities,  which  do  not  necessarily  land  on  the 
desired  paths,  as  the  most  probable  reason  for  the  failure  of  schemes  to  converge  to  the 
correct  weak  solution  on  the  desired  paths,  the  proposed  subcell  resolution  technique  within 
a  Runge-Kutta  time  discretization  seems  to  have  robustness  problems  for  very  strong  dis¬ 
continuities,  especially  during  interactions  of  such  discontinuities.  An  improvement  on  the 
robustness  of  the  algorithms  when  transition  points  are  removed  or  significantly  reduced 
constitutes  our  future  work. 

Our  discussion  is  restricted  to  the  one-dimensional  case.  Preliminary  study  on  the  ex- 
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Figure  8.10:  Density,  velocity  and  pressure  for  Example  10.  t  =  0.0001.  Solid  line:  the  exact 
solution.  Symbol:  the  numerical  solution. 
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tension  of  our  methodology  to  two-dimensional  problems  (not  reported  in  this  paper),  im¬ 
plemented  in  a  dimension  by  dimension  fashion  (not  dimension  splitting),  has  indicated 
limitations  for  some  of  the  complex  discontinuity  wave  patterns.  This  may  be  due  to  the 
complication  of  two-dimensional  waves  and  the  limitation  of  the  subcell  resolution  technique 
to  capture  accurately  the  correct  states  and  the  corresponding  paths  at  the  discontinuities. 
Further  investigation  of  effective  multi- dimensional  algorithm  in  this  context  also  constitutes 
our  future  work. 
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